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We present a new dynamic off-equilibrium method for the study of continuous transitions, which 
represents a dynamic generalization of the usual equilibrium cumulant method. Its main advantage 
is that critical parameters are derived from numerical data obtained much before equilibrium has 
been attained. Therefore, the method is particularly useful for systems with long equilibration times, 
like spin glasses. We apply it to the three-dimensional Ising spin-glass model, obtaining accurate 
estimates of the critical exponents and of the critical temperature with a limited computational 
effort. 
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Monte Carlo simulations combined with finite-size scal¬ 
ing (FSS) methods are, at present, the most success¬ 
ful tool for the identification of continuous transitions 
and the determination of the critical parameters. In 
this approach there are two main obstacles to a pre¬ 
cise determination of the critical parameters. On one 
side, scaling corrections, related to the subleading irrel¬ 
evant renormalization-group (RG) operators, may mask 
the asymptotic critical behavior, which shows up only 
when the system size L becomes large. On the other side, 
the Monte Carlo dynamics becomes increasingly slow as 
the critical point is approached, so that thermalization 
and equilibrium autocorrelation times become large as L 
increases, hampering large-L simulations. These prob¬ 
lems are particularly serious in systems with quenched 
disorder. They occur, for instance, when studying the 
paramagnetic-glassy transition in the ±J Ising model 
with random ferromagnetic and antiferromagnetic cou¬ 
plings, which represents a standard theoretical labora¬ 
tory to understand the effects of quenched disorder and 
frustration. Its Hamiltonian is given by [Ij 

H = ^ ^ Jxy^x^yi ( 1 ) 

Gv) 

where the sum runs over all nearest neighbors {xy) of 
a cubic lattice, = ±1 are Ising spins, and Jxy are 
quenched random couplings that take the values ±1 with 
equal probability. At the transition and, even worse, in 
the low-temperature phase, the standard Metropolis dy¬ 
namics is very slow, so that equilibration becomes very 
difficult, even for relatively small systems. Moreover, 
equilibration times depend strongly on the disorder re¬ 
alization, so that a sample-by-sample analysis is needed 
to guarantee that all measurements are obtained from 
well-equilibrated sarnples 0. At present, even by using 
dedicated machines [3|, it seems impossible to go much 
beyond sizes L = 30-40. 

In this work we discuss a dynamic method to deter¬ 
mine the critical temperature and the critical exponents. 


We will discuss it in the context of the ± J Ising model, 
but the method and the results are completely general, 
so it can be applied to the study of any continuous tran¬ 
sition in pure or disordered systems. The method, which 
does not require the system to be in equilibrium, has 
two advantages. First, the difficult and time consuming 
(at least in disordered systems) task of verifying equili¬ 
bration is no longer needed. Second, we can stop the 
simulation much before equilibrium has been reached, 
saving a considerable amount of computing time. The 
method we discuss is somewhat different from previous 
off-equilibrium methods (see, e.g.. Refs. 0-12 and refer¬ 
ences therein). Indeed, in most of those works it is gen¬ 
erally assumed that L is so large that finite-size effects 
are negligible, a condition that is easily satisfied in pure 
systems but not in the random case. On the contrary, 
we will use the finite-size dependence of physical observ¬ 
ables to estimate the critical parameters. Our method is 
essentially an off-equilibrium generalization of the usual 
Binder crossing method. 

As in FSS equilibrium computations, we begin by con¬ 
sidering a RG invariant quantity C/(t,L,/3) as a function 
of the Monte Garlo time t, inverse temperature /3, and 
system size L. According to RG, for L and t large and 
close to the critical point /3c, scales as 13-1^ 


U{t, L, P) — f}i{tL ^, e) -\- Ui^(/3)L ^,e)-|-... 

( 2 ) 

where next-to-leading scaling corrections have been ne¬ 
glected. Here w is the leading correction-to-scaling expo¬ 
nent, ittc(/3) the associated nonlinear scaling field satisfy¬ 
ing u^(/3c) ^0, e = , where rv /3 - /3c is the 

temperature nonlinear scaling field which parametrizes 
the distance from the critical point. Equation ([2]) de¬ 
pends also on the dynamic critical exponent z, which 
represents an additional parameter to be determined. To 
avoid any reference to z, we now reparametrize the time 
evolution in terms of the correlation length or, better, 
in terms of the RG invariant ratio = £,!L. When us- 
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ing an initial disordered configuration, ^ is an increasing 
function of t, any function of tL~^ can be equivalently 
reexpressed in terms of R^, so that we write 

U{t, L, P) = fuiR^, e) + u^{P)L-^gu{Ri.e) + ..., (3) 

which is defined for < i?^^eq(e), where i?^_eq(e) is the 
equilibrium value of R^ for the given e. Equation ([S]) is 
the basic relation we use to compute critical tempera¬ 
ture and exponents. Indeed, ignoring scaling corrections, 
close to the critical point Eq. ([3|) can be expanded in e, 
obtaining 

U{t, L, /3) = fu{Ri, 0) + (/3 - 0) +... (4) 

At fixed R^, the quantity U{t,L,P) behaves exactly as 
in the equilibrium case: /3c is determined as the crossing 
point and v is obtained by computing the slope at Pc- 
However, in this formulation equilibration is not needed. 
Equation (|31) is valid for any value of R^, so one might 
think of choosing a small value for such a parameter, 
reducing significantly the length of the runs. However, 
one must not forget that the method is intrinsically a 
finite-size method; hence, it can only work if finite-size 
effects are not too tiny, and this, in turn, requires R^ to 
be not too small. Mathematically, these considerations 
can be understood by considering Eq. ®. The method 
is expected to be precise if the coefficient 0) is not 

too small. Such a coefficient depends on R^ and it is ex¬ 
pected to increase with R^. In particular, it is expected 
to be very small for R^ small, so that, if one chooses a 
small value of R^, the crossing becomes undetectable, un¬ 
less statistical errors are very tiny. Hence, R^ should be 
chosen small, but still large enough to have a reasonable 
sensitivity of the results on system sizes. 

To validate the method, we apply it to the determina¬ 
tion of the critical point and of the critical exponents in 
the ± J Ising model. We perform large-scale simulations 
on cubic lattices of volume L^, with 8 < L < 64, consid¬ 
ering five values of P between 0.880 and 0.910. Statistics 
is a crucial factor in the analysis and hence we consider 
a very large number Ng of samples for each L and /3. 
Typically, Ng varies between 3 • 10® and a few million. 
Only for L = 48, 64 is Ng smaller: Ng = 6 ■ lO'^, 10"^ in 
these cases. Essentially all runs end when the system is 
still out of equilibrium. In most of the cases, data extend 
only up R^ ~ 0.5, in some cases even less (at equilibrium 
Q = 0.652(3) at the critical point). Simulations were 
performed on a small GPU cluster using a very efficient 
asynchronous multispin coding technique [l^ . In each 
run we simulate together 32fc different disorder realiza¬ 
tions with four replicas for each disorder realization. The 
value of k is tuned for each L to have the best perfor¬ 
mance of the GPUs. As a result, one spin flip takes 2.9 
ps (essentially for all sizes) on the GTX Titan, the fastest 
GPU we have. The simulations presented here took ap¬ 
proximately 3.1 CPU years of the GTX Titan GPU. 


To apply Eq. m we must define the quantities U. We 
consider 5 different Binder cumulants defined in terms of 
the overlap between two different replicas Qx = 
Moreover, we must somehow parametrize the scaling 
functions fu{R^, e) and guiR^, e)- Since the data belong 
to a small temperature interval, we use the expansion (|3]) 
to first order in e. We have also performed some analyses 
using a second-order approximation, without observing 
significant differences. As for the correction-to-scaling 
function, we have verified that we can assume it to be 
independent of temperature. Finally, we should make 
approximations for the nonlinear scaling fields. Relying 
on the analysis of Refs. 00, we set U/3(/3) = P — Pc 
and Ucj{P) = 1, neglecting the additional corrections. 
Given the small temperature interval we consider, these 
approximations should hold quite precisely. Hence, each 
U{t, L, P) was fitted to 

C7(t, L, P) = Pi{R^) + P2(Rc)(/ 3 - /3,)Li/" + P:i{R^)L-'^, 

( 5 ) 

with Pi(i?^), Pi(P^), and Pz{R{) polynomials of degree 
6, 3, and 3, respectively. The fit of the five renormalized 
couplings is quite complex, as we take w. Pc, v, and the 
coefficients of the polynomials as free parameters. As a 
whole, there are 78 free parameters that must be opti¬ 
mized. In the fits we have not taken into account the time 
correlations among data at the same P and L, so that sta¬ 
tistical errors (computed using the jackknife method) are 
not a priori optimal. To verify that such neglect is not 
relevant for the final estimates, we have performed some 
fits of a single cumulant taking time correlations into ac¬ 
count. The corresponding estimates and error bars are 
essentially equal to those obtained without including sta¬ 
tistical correlations. 

As usual in this type of analyses, the most difficult 
issue is the estimation of the systematic errors due to 
the neglected correction terms. This is very important 
here, since the attainable values of L are quite small. 
We have thus performed fits with several types of cuts. 
We perform fits including each time only data satisfy¬ 
ing L > Pmin, c > ^min, and R^ > Considering 

several values for Pmin, Cmin, and Results ob¬ 

tained by taking 3 < ^min < 5, 8 < Pmin < 12, and 
0 < P{,min < 0.4 show some scatter, which is some¬ 
what larger than statistical errors, indicating that the 
neglected systematic effects may be as important as the 
statistical ones. The most crucial parameter is ^min- 
When such a parameter is increased from 3 to 4, the ex¬ 
ponent w decreases sharply, by more than one error bar, 
while Pc increases. Such a systematic drift occurs also 
when ^min is further increased to 5, but now the change 
is much less than one error bar. Therefore, the results 
we quote correspond to fits with ^min = 4. For such a 
value of ^min we obtain Pc = 0.911(2), 0.916(4), 0.909(4) 
for Pmin = 8,10,12 and P^.min = 0, and Pc = 0.911(2), 
0.909(2), 0.909(3) for P^^min = 0, 0.2, 0.4 and Pmin = 8. 
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FIG. 1: Plot of AC /4 versus t = {P - for = 0.25 

(top), = 0.35 (middle) and = 0.45 (bottom). We set 
Pc. = 0.910, uj = 1.3 and v = 2.47. Here U 4 is the standard 
Binder cumulant, as defined in Refs. 1.0. Symbols: empty 
square {L = 8 ), empty circles {L — 10), empty triangles {L — 
12), empty pentagons (L = 16), filled squares (L = 20), filled 
circles {L = 24), filled triangle {L = 32), filled pentagon (L = 
48). 


No systematic trends can be observed, all estimates be¬ 
ing consistent within errors. Except for one estimate, 
all results (with their errors) we are quoting here lie in 
the interval 0.906 < Pc ^ 0.913. Therefore, we take 
Pc = 0.910(4) as our final estimate. The error, which 
is twice the error affecting the results with Tmin = 8, is 
somewhat subjective and should take into account the ef¬ 
fect of the neglected next-to-leading scaling corrections. 
Analogously, we can estimate w and v obtaining 

a; = 1.3(2), :/ = 2.47(10). (6) 

The estimates of cu are strongly correlated with those of 
Pc- the larger pc, the smaller w is. If /3c = 0.906, fits 
keeping Pc fixed give uj ~ 1.5, while w « 1.1 is obtained 
by fixing Pc = 0.914. The exponent v is instead much 
less correlated with Pc, changing at most by ±0.03 when 
Pc varies by ±0.004. 

To show the quality of the results, in Fig. [I] we report 
AC/ 4 , defined by 

AC/4(/3, L, R^) = UiiP, L, R^) - P3{R^)L-‘^, (7) 

versus e. We consider R^ = 0.25, 0.35, and 0.45. Very 
good scaling is observed, confirming the correctess of the 
scaling Ansatz and the accuracy of the estimates of the 
critical exponents. Note also that the data lie on an 
essentially straight line, validating our choice of expand¬ 
ing fu(R^,e) to first order in e. From the figure, we 
can also clarify why a large number of samples, of or¬ 
der 10®, is needed to estimate the critical parameters. 
For instance, C /4 at R^ = 0.35 varies by 0.04 within our 
temperature interval. Therefore, the temperature depen¬ 
dence of the data can be observed only if the errors on C /4 


are significantly less than 10“^, for instance, if they are 
equal to 10“®. Since errors scale as a/y/Wl with a ~ 1 
for all values of L, a 10“® error is obtained by taking 
Ng ~ 10®. Note that this requirement is not specific 
of the off-equilibrium method we use. Also equilibrium 
analyses require Ng to be large a 0,1111. 

It is interesting to compare these results with previ¬ 
ous ones, see Table H] (older estimates are summarized 
in Ref. l23h . For the critical-point position, our estimate 
Tc = 1/Pc = 1.099(5) agrees within errors with the esti¬ 
mates Tc = 1.102(3) and Tc = 1.109(10) of Refs. aim, 
obtained from the analysis of equilibrium results. Our 
error is larger than that reported in Ref. a but note that 
our final error includes a subjective estimate of the sys¬ 
tematic error. Had we reported only the statistical error 
for Lmin = 8, we would have obtained the same accuracy. 
The estimates of v are also consistent, while our final esti¬ 
mate of w is slightly larger, though still consistent within 
error bars, th an p revious ones. The off-equilibrium esti¬ 
mates of Ref. [ill are consistent with ours, but less pre¬ 
cise. Previous dynamic estimates of Tc are instead not 
consistent. It is now clear that the reported errors are 
underestimated, as a consequence of the neglect of the 
subleading scaling corrections in the analyses. 

The method we have discussed represents a significant 
improvement with respect to equilibrium analyses. In¬ 
deed, since the scaling variable is tL~^, the time needed 
to extend Metropolis runs from any value of R^ to equi¬ 
librium scales as L^, i.e., as given that z Ri 7 for the 
Ising spin glass. Therefore, the advantage is very large 
and increases rapidly with L. To make a fair comparison 
with equilibrium studies, we should, however, take into 
account that in those studies one combines the parallel¬ 
tempering method 2^ with the Metropolis or heat-bath 
algorithm. It is not clear how equilibration times scale 
for this combined algorithm, and in particular, how long 
it takes to thermalize the hard samples. However, the 
results reported in Ref. |3 are consistent with a sample- 
dependent time that scales as for the samples that 
equilibrate fast and as for those that are slower. The 
off-equilibrium method is still significantly faster. A more 
direct comparison can be obtained by using the results 
published in Ref. 0. In our simulations at the critical 
point, runs extending up to R^ ~ 0.5 require 2.5 • 10®, 
16 • 10® Metropolis sweeps for L = 24 and 32, respec¬ 
tively. In the parallel-tempering simulations for L = 32 
of Ref. 1^ the number of iterations discarded for thermal- 
ization varies between 8 • 10® and 500 • 10® (on average 
13 • 10®) sweeps. Taking into account that 22 systems 
at different temperature are simulated together, our sim¬ 
ulations are shorter by a factor of 10 at least. If one 
were stopping the off-equilibrium runs at R^ = 0.40, one 
would gain an additional factor of 3 for this value of L. 

In spite of the significant improvement with respect 
to equilibrium studies, the computing time needed for 
a simulation scales as even off-equilibrium, since we 
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TABLE I: Estimates of Tc and of the critical exponents by off-equilibrium methods. Results of Refs. are obtained in 

equilibrium simulations. The exponents /3 and 7 are related to v and 77 by /3 = u(l -|- ri)/2, 7 = v(2 — rf). 



P 

V 

UJ 

z 

p 

7 

zv 


Ref. 4 

1 . 12 ( 12 ) 



5.7(5) 

Ref. ^ 

1.17(4) 

1.5(3) 


6.2(2) 3.6(6) 

Ref. ^ 

1.19(1) 


- 0 . 22 ( 2 ) 

5.7(2) 

Ref. 7 

1.154(3) 



0.52(9) 

Ref.^ 




6.86(16) 

Ref. ^ 

1.135(5) 




Ref. JT 

1.18(1) 

1.40(5) 

- 0 . 20 ( 1 ) 


Ref. JT , 22 

1.115(15) 

2.2(3) 

-0.380(7) 

6.79(6) 

Ref. 12 




5.85(9) 

this work 

1.099(5) 

2.47(10) 

-0.39(1) 1.3(2) 

6.80(15) 

Ref. 21 

1.109(10) 

2.45(15) 

-0.375(10) 1.0(1) 


Ref. 3 

1.1019(29) 2.562(42) 

-0.3900(36) 1.12(10) 




tL-^ 


FIG. 2: Plot of C/4 and versus tL ^ for P = 0.910. We set 
a = 6.80. 


need to collect data at fixed for all values of L. This 
requirement makes our method not suitable to investigate 
large system sizes. Since errors are independent of system 
size as a consequence of the absence of self-averaging, the 
Monte Carlo time needed to obtain the same statistical 
errors scales also as L^. This explains why we have not 
considered lattices with T > 64. If we increase L, we 
should increase t at the same time, making simulations 
far too long. 

The analysis we have performed for the renormalized 
couplings can be extended to the susceptibility. The 
finite-time scaling behavior can now be written as 

Inx = {2-77)lnL + Pi(i?5) + (/3-^,)Li/"P2(i?«) + 

L-^P3{Ri) + Pm. ( 8 ) 

where the last term Pi{fi) is the contribution of the non¬ 
linear scaling field associated with the magnetic field, see 


Ref. 0, [U for a discussion. A good parametrization 
is obtained by taking Pi{R^), P 2 {R^), P 3 {R^) as poly¬ 
nomials of degree 6 , 3, 3, respectively, as before. For 
the P 4 {I 3 ), we set P 4 = 04 /?. We obtain the final esti¬ 
mate T] = —0.39(1), which is fully consistent with those 
of Refs. sm 


Finally, we estimate z by requiring data to satisfy the 
general scaling form ©■ We obtain 2 ; = 6.80(15), where 
the error should be quite conservative. The scaling be¬ 
havior of R^ and U 4 is shown in Fig. [5J Scaling correc¬ 
tions are clearly visible, but large L data appear to fall 
onto a single universal curve as L increases. The value we 
obtain is in agreement with the estimate of Refs.Js, 11 at 
T = 1.1. Instead, it is larger than those of Refs. [3 [alS. 
However, note that in all these works no scaling correc¬ 
tions, crucial to control possible systematic errors, were 
included in the analyses (they play a fundamental role in 
the derivation of our result). 


Let us now summarize our results. We have presented 
a new dynamic off-equilibrium method suitable for the 
determination of the critical temperature and of the crit¬ 
ical exponents. Such a method represents a significant 
improvement with respect to previous ones. In particu¬ 
lar, there is no need for L to be large enough to avoid 
finite-size effects—thus, a source of systematic errors is 
absent—nor does it require an a priori knowledge of the 
critical temperature. We have used the method to deter¬ 
mine critical exponents and temperature for the ± J Ising 
model. With a relatively modest investment of comput¬ 
ing time, thanks also to a very efficient GPU multispin 
code, we obtain results that have a comparable preci¬ 
sion with that of the estimates of Ref. Q , which are the 
most precise equilibrium estimates available today. The 
method is completely general and can be applied to any 
pure or disordered system. 
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